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Abstract 



We use a phenomenological field theory, reflecting the symmetries and con- 
servation laws of sandpiles, to compare the driven dissipative sandpile, widely 
studied in the context of self-organized criticality, with the corresponding 
fixed-energy model. The latter displays an absorbing-state phase transition 
with upper critical dimension d c = 4. We show that the driven model ex- 
hibits a fundamentally different approach to the critical point, and compute 
a subset of critical exponents. We present numerical simulations in support 

of our theoretical predictions. 

PACS numbers: 64.60.Lx, 05.40. +j, 05.70.Ln 



1 



Typeset using REVTpjX 



A wide variety of nonequilibrium systems display transitions between "active" and "ab- 
sorbing" states: examples are epidemic processes |1| catalysis 0, directed percolation (DP) 
||, and the depinning of interfaces in quenched disorder |3J]. When driven continuously, such 
systems may exhibit stick-slip instabilities, or broadly distributed avalanches, commonly as- 
sociated with self-organized criticality (SOC) |5|J§]. 

SOC sandpiles f| possess an infinite number of absorbing configurations (i.e., from which 
the system cannot escape), and are placed, by definition, at the critical point in a two- 
dimensional parameter space 0H resembling that of directed percolation (DP) || or contact 
processes P-pH- 

Under an external drive (i.e., input of particles at rate h), the system jumps among ab- 
sorbing configurations via avalanche- like rearrangements. Close to the absorbing-state phase 
transition, a slow drive induces avalanches whose size distribution decays as a power law - 
the hallmark of SOC. What distinguishes the sandpile from other models with absorbing 
states is a conservation law: avalanche dynamics conserves the number of grains of "sand," 
and the order parameter is coupled to this conserved field ||. 

In this Letter we use a phenomenological field theory of sandpiles to show how conser- 
vation alters the phase transition. The critical behavior for h — > (the SOC limit) differs 
from that for h = 0. In particular, when driving and dissipation are absent, the sandpile 
shows an absorbing-state phase transition (with d c — 4). Our approach clarifies the effect 
of driving on dynamic phase transitions, and resolves several long-standing issues regard- 
ing sandpiles, such as the upper critical dimension, the effect of conservation on critical 



exponents, and universality classes fll2|-|16|j. We perform extensive simulations to check our 
theoretical predictions. 

In sandpiles ||, each site i of a d— dimensional lattice bears an integer variable Zi > 0, 
which we call energy. When a site reaches or exceeds a threshold z c it topples: — > Zj — z c , 
and Zj —>■ Zj + 1 at each of the g nearest neighbors of i. Energy is fed into the system at rate h, 
and is dissipated at rate e during toppling [[?[]. At each time step, each site has a probability 
oc h to receive an energy grain; in each toppling, a grain is lost with probability oc e. These 



rules generalize the original Bak, Tang and Wiesenfeld (BTW) sandpile automaton ||, which 
is recovered in the limit h — > and e — > ]7J§|. While the BTW model restricts dissipation 
to the boundaries, we focus on (conceptually simpler) bulk dissipation; most conclusions 
apply to the boundary dissipation well. We also consider the Manna sandpile, in 

which z c = 2 and two neighboring sites are chosen at random to receive energy | |17| . 

In the slow-driving limit, each energy addition is followed by an avalanche of s topplings; 
the avalanche distribution has the scaling form P(s) = s~ T G(s/s c ), where the cutoff scales 
as s c ~ £ D . The correlation length £ scales with dissipation as £ ~ e~", and is related to the 
characteristic avalanche duration t c ~ £ z . 

The order parameter is p a , the density of active sites (i.e., whose height z > z c ) f7|||; if 
p a = the system has reached an absorbing configuration. In a coarse-grained description, 
we study the dynamics of a local order-parameter field p (x, t), bearing in mind that the 
energy density C{x,t) is (f° r e = /i = 0), a conserved field. Variations of the local energy 
density are due to: (i) the external field, h; (ii) dissipation attending toppling: — ep a ; (iii) 
a diffusion- like contribution: (1 — e)V 2 p a . The latter arises because a gradient in p a leads 
to a current: the excess in the mean number of particles arriving at x from the left, over 
those arriving from the right, is j x (x,t) = —(1 — e)d x p a . The net inflow of particles at x is 
therefore —V • j = (1 — e)V 2 p a - Defining an energy diffusion constant oc 1 — e, we write 
the continuity equation for the energy density, 

dC(x t) 

y ^ = A:V 2 Pa(x, t) - ep a (x, t) + h(x, t) + % (x, t), (1) 

where the driving field h(x, t) = h + r)h(x, t), with h a nonfluctuating term and 7]^ zero-mean, 
uncorrelated Gaussian noise. The last term is dynamically generated Reggeon-field-theory- 
like (RFT) noise |18] f^(x, t) ~ y / p fl (x, t)?y(x, t), with rj uncorrelated Gaussian noise. This 
term vanishes, as it must, in the absorbing state, p a = 0. 

The equation for the order-parameter field is readily obtained || by extending the mean- 
field theory (MFT) of Ref. (7J. With p e (x,t) the local density of "critical" sites (i.e., with 
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height z c — 1), we have 

dp a (-X,t) 



D a V 2 pa(x,t) - p a (x,t) 



Of 

+ (9- e)Pa(x, t)p c (x, t) + /ip c (x, t) + 7? a (x, t), (2) 



where r/ a (x, t) is a RFT-like noise whose amplitude is proportional to Jp a (pc,t). The first 
two terms represent toppling jl9 |; the terms oc p c represent critical sites becoming active 



upon receiving energy, whether from the external drive, or from toppling neighbors. 

In the stationary state, we can avoid ensnarement in an infinite hierarchy of equations 
||, by eliminating p c in favor of ( and p a . In the Manna model, we simply invoke nor- 
malization: p c = C — Zapa, where z a > 2 is the mean height of active sites. For BTW we 
use the phenomenological Ansatz: p c (x, t) = [1 — p a (x, £)]/[£(x, t)]. That is, the fraction / 
of nonactive sites that are critical can be expressed as a single-valued function of the local 
energy density. In the slowly-driven stationary state, ( ~ ( s and / ~ p£°, where ( s and 
p^° are the stationary average values of the energy and the critical-site density. We expand 

/(£) = pf + AA( H , where A( = C(x, t) - ( s , and A > 0. We test the validity of these 

assumptions by simulating the two-dimensional BTW model on a lattice of 80 x 80 sites, 
at ( = ( s = 2.125. To determine /(C), we measure the average energy, and active- and 
critical-site densities in cells of 10 x 10 sites. The conditional probability density P(/|C) is 
unimodal and peaked (see Fig.l); the mean increases linearly with (, indicating that p c is 
well-approximated by (1 — p a )f{(), with / linear in the neighborhood of ( s . 

Eqs. (P and @, describing the coarse-grained dynamics of sandpiles, resemble the field- 
theory for the pair contact process (PCP) ||11|| , another model with infinitely many absorbing 
states. As in the PCP, when h = all configurations £(x) consistent with p a = are 
absorbing, and the order parameter is coupled to a non-order-parameter field playing the 
role of an effective creation rate. The essential difference in the sandpile is that the field ( 
is conserved. In the following we consider separately the cases of slow driving (h — > + ), 
corresponding to the SOC sandpile, and of fixed energy: h = and e = 0. 

(i) Driven sandpile: The system attains its stationary state by the very slow addition 



of energy. In this limit [h — > + ), a complete time scale separation, between toppling, on 
one hand, and addition and dissipation, on the other, sets in |7|||,[J. In the stationary state, 
energy balance forces a subset of the critical exponents to take their mean-field (MF) values 
in any spatial dimension d, as we now show. Integrating Eq. ([!]) over space and averaging 
over the noise yields (p a ) = V~~ l J d d x{p a (]t,t)) = h/e 0, which implies that the zero-field 
susceptibility \ = dp a /dh, diverges as e~ 7 with 7 = 1. Taking the functional derivative of 
Eq. (fj) with respect to h(pd), and averaging over the noise, we obtain an equation for the 
static response function Xe( x — x ') — (^Pa(x)/5/i(x')): 

- D c V 2 Xe (x - x') + e Xe (x - x') = (5(x - x'), (3) 

which yields, for large r, Xe( r ) oc r 2 ~ d e~ r ^, where the correlation length £ ~ e~ u , with 
v = 1/2. These results depend solely upon stationarity and translation invariance p0| . 



Although the exponent values coincide with the MF ones, they have not been obtained by 
MF arguments, and are valid beyond MFT, as confirmed by simulations in 2 < d < 6 fT4|| . 

While the remaining exponents are in principle also determined by Eqs. ([I]) and (|2|), a 
full analysis, involving the double limit t — > 00, h — ► + (the order cannot be interchanged), 
promises to be a knotty problem. The critical properties emerge as e — > 0, which must be 
taken subsequent to the above limits, since a stationary state demands e > h. The upper 
critical dimension, however, can be found by power-counting analysis. The evolution of a 
localized perturbation p a around the slowly-driven stationary state (h/e = 0) is given by 

do fx t) 

K m J = ^aV 2 pa(x, t) - rp a (x, t) + pAC(x, t)p a (x, t) - M p 2 (x, t) + Va {x., t), (4) 

where r ~ e and p and u are coupling terms generated by the elimination of p c in favor of 
(. We can consider in Eq.s (|l|)and (JD the usual rescaling x — ► bx', t — > b z t' and p a — > b 5a p' a , 
and the rescaling of the energy field £ — ► fe^C'. The rescaled coupling constants show a 
MF fixed point for r = e = and z — 2. In this case the nonlinear terms couplings scale 
as p ~ u ~ 6 4_d . Thus, nonlinear terms are irrelevant when d > 4, defining an upper 



critical dimension d c = 4 pl| . Since r = e = is the fixed point, the dissipation rate is 



the (temperature-like) control parameter, with critical value e c = 0, emphasizing the role 
of conservation in slowly-driven sandpiles. By studying Eq. (^) in the slowly-driven MF 
stationary state, (A£ = p a = 0, pf = 1/g, and neglecting noise terms) we also obtain the 
MFT exponents for avalanche spreading: D = 4, r = 3/2 and z = 2, in agreement with 
earlier analysis J7J. 

(ii) Fixed- energy sandpile (FES): When h = e = 0, the total energy / o((x, t) is 
conserved and plays the role of a control parameter. In this case, Eq. (|]) reduces to 
d£(x,t)/dt = V 2 p a (x, t) + ?7f (x, t), where 7fc(x, £) is a conserved noise. Substituting the 
formal solution of this equation into Eq. (§) yields 



dt 



D a V 2 pa(x, t) - r(x)p a (x, t) - 6(xK(x, t) 



wp a (x,t) [ V 2 p a (x,t')dt' + rj a (x.,t), (5) 

JO 



where we neglect higher-order terms, irrelevant by naive power- counting analysis. The 
coefficients r and b depend on position through on the initial value of C( x )- Eq. ([5]) is 
the Langevin equation of RFT, save the non-Markovian term and the spatial variation of 
r and b; d c = 4, as in RFT. In MFT, replacing £(x, 0) by the spatially uniform (, we have 
dpajdt = —rp a — bp 2 al i.e., the MFT of directed percolation (DP), with critical point at r = 0, 
fixing, in turn, the critical energy density, ( c . Close to the critical point (|A£|/£ C << 1), 
r ~ A(. For d < d c , the critical fixed point will be renormalized to r = r*, defining a 
renormalized ( c . Above ( c , we have an active stationary state with p a ~ (A£)^; for ( < ( c , 
the system falls into an absorbing configuration in which p a = 0. 

Thus the FES approach to criticality (( — > ( c ) is fundamentally different from the driven 
case, (h,e — > 0, followed by h/e — > 0). Note that ^ is lost as an independent parameter 
once h and e are nonzero (Slow driving pins ( at its critical value: if it exceeds ( c , activity 
is generated, and thereby dissipation). The behavior at the critical point is described by 
our theory with h — e — and £ = £ c . As in other models with infinitely many absorbing 
configurations, the avalanche behavior depends intimately on the initial configuration. It is 
also worth remarking that in the stationary driven case, the dynamics can explore only a 
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set of recurrent configurations ||22|| . The FES may instead explore transient configurations 
that could account for the different critical behavior. 

To better understand its scaling, we simulated the FES with statistically homogeneous 
initial conditions We considered the BTW model with periodic boundary conditions at 
(e = 0,h = 0). Initial configurations are generated by distributing at random a fixed 
number iV of particles among the L d sites. This defines the strictly conserved energy density 
( = N/L d . Once all N particles have been placed, active sites topple at a unit rate with a 
sequential updating rule. We studied the transition from the active to the absorbing states 
as we varied (. In d = 2, using system sizes extending up to L = 1280, we find ( c = 2.125, 
(3 = 0.59(1), v = 0.79(4), and z = 1.74(4) p3fl . (Figures in parentheses denote uncertainties.) 
The corresponding DP exponents are 0.583(4), 0.73(2), and 1.76(3). Simulations of the four- 
dimensional model yield ( c = 4.11(1) and (3 = 1.00(1), in good agreement with theoretical 
results, which predict MF values in d > d c = 4 (see also Fig.2). 

These results are compatible with the DP universality class, suggesting that the non- 
Markovian term is irrelevant, at least for homogeneous initial conditions. Preliminary results 
of direct integration of Eq.(^) indicate DP-compatible behavior for homogeneous initial 
conditions; the non-Markovian term does appear to alter the spreading exponents, as in 



other multiple-absorbing state models |3[TT]]. O n the other hand, we find analytically that 



the non-Markovian term is relevant at the RFT fixed point below d = 4, and it has to be 
taken into account to determine the asymptotic scaling properties. This implies that to fully 
understand the effect of this term we need to perform a full RG perturbative expansion as 
well as larger numerical simulations. 

In summary, our field theory elucidates the effect of driving on the critical behavior of 



sandpiles through their connection with other absorbing-state phase transitions |24| . Field- 
theoretic analysis shows that driving does not cause such a change of critical behavior in 
the PCP, as we have verified numerically. We believe the crucial difference is the absence 
of a conservation law in the PCP. Indeed, a different kind of conservation law ( "local parity 
conservation" ) also changes the critical behavior of contact-process-like models |25| . 



Finally, we remark that our framework applies equally to the BTW and Manna sandpiles, 
even though the latter has a stochastic toppling rule [|17[]; the additional noise generates no 
further relevant terms. Though it remains an open question in the context of simulations [f26[ . 
recent large-scale numerical studies and rigorous arguments support the shared universality 
of the two models |27 . 

Our approach suggests several paths for further investigation. An open question concerns 
the critical behavior of fixed-energy models with nonhomogeneous initial conditions. A 
full renormalization group treatment of the field theory, while challenging, should yield 
systematic predictions for avalanche exponents. Analysis of different boundary conditions 
should lead to a better understanding of the scaling anomalies exhibited by sandpile models. 

During the completion of this work, we learned of large scale simulations of the BTW- 
FES model with parallel updating performed by P. Grassberger. The results indicate that 
deviations from DP behavior persist at large system sizes, suggesting that the non-Markovian 
operator in our field-theory could become relevant below d c . 
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ments and discussions. M.A.M., A.V. and S.Z. acknowledge partial support from the Eu- 
ropean Network Contract ERBFMRXCT980183. M.A.M. is partially supported from the 
TMR European Program ERBFMBIGT960925. 
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FIG. 1. Fraction of critical nonactive sites, / (shifted by p c ) versus C~Cc for the BTW in d = 2. 
Inset: conditional probability density at £ = 2.12 for the same case. 
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FIG. 2. Stationary active-site density as a function of £ — £ c for the fixed-energy BTW model 
in d = 4. The inset shows p a in a larger range of £ — £ c and for different lattice sizes L. 
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